Importing and
organizing normalized IMR90 gene expression data
IMR90_edat<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/imr90_wgcna/IMR90_edat.RDS")
Taking mean
expression value per condition per time point
IMR90_ER<-IMR90_edat[, grep("ER", colnames(IMR90_edat))]
IMR90_GFP<-IMR90_edat[, grep("GFP", colnames(IMR90_edat))]
ER_time_point<-c("ER_0", "ER_4", "ER_8",
"ER_12", "ER_16", "ER_20",
"ER_24", "ER_32", "ER_40",
"ER_48")
GFP_time_point<-c("GFP_0", "GFP_4", "GFP_8",
"GFP_12", "GFP_16", "GFP_20",
"GFP_24", "GFP_32", "GFP_40",
"GFP_48")
IMR90_ER_mean<-data.frame(genes = rownames(IMR90_ER))
for (i in 1:length(ER_time_point)){
time_point<-ER_time_point[i]
sub_df<-as.data.frame(IMR90_ER[, grep(time_point, colnames(IMR90_ER))])
sub_df[,"mean"]<-apply(sub_df, 1, mean)
df_mean<-as.data.frame(sub_df[, "mean"])
colnames(df_mean)<-time_point
IMR90_ER_mean<-cbind(IMR90_ER_mean, df_mean)
}
IMR90_GFP_mean<-data.frame(genes = rownames(IMR90_GFP))
for (i in 1:length(GFP_time_point)){
time_point<-GFP_time_point[i]
sub_df<-as.data.frame(IMR90_GFP[, grep(time_point, colnames(IMR90_GFP))])
sub_df[,"mean"]<-apply(sub_df, 1, mean)
df_mean<-as.data.frame(sub_df[, "mean"])
colnames(df_mean)<-time_point
IMR90_GFP_mean<-cbind(IMR90_GFP_mean, df_mean)
}
IMR90_ER_mean_matrix<-as.matrix(IMR90_ER_mean[, 2:11])
IMR90_GFP_mean_matrix<-as.matrix(IMR90_GFP_mean[, 2:11])
IMR90_ER_GFP_DELTA<-IMR90_ER_mean_matrix-IMR90_GFP_mean_matrix
IMR90_ER_GFP_DELTA<-as.data.frame(IMR90_ER_GFP_DELTA)
rownames(IMR90_ER_GFP_DELTA)<-IMR90_ER_mean$genes
Adjusting gene
expression value of ER samples by substracting the value of GFP samples
at the same time point
IMR90_ER_GFP_DELTA_df<-reshape2::melt(as.matrix(IMR90_ER_GFP_DELTA))
colnames(IMR90_ER_GFP_DELTA_df)<-c("genes", "time", "ER-GFP")
IMR90_ER_GFP_DELTA_df$time<-unlist(lapply(as.character(IMR90_ER_GFP_DELTA_df$time), function(x) strsplit(x, "_")[[1]][2]))
IMR90_ER_GFP_DELTA_df$time<-as.numeric(IMR90_ER_GFP_DELTA_df$time)
IMR90_ER_GFP_DELTA_df$`ER-GFP`<-as.numeric(IMR90_ER_GFP_DELTA_df$`ER-GFP`)
Plotting the gene
markers for MCC and neuroendocrine tumors in ER 48 hours samples
# gene markers
neuroendocrine_markers<-c("ENO2", "NEFM", "NEFH", "NMB", "HES6")
neural_genes<-c("EFNB1","SEMA4D", "SLIT2")
keratin_markers<-c("KRT8", "KRT18")
wnt_down_gene_list<-c("WNT5A","WNT5B", "FZD2", "FZD7","FZD8", "WNT16", "TCF7L2")
wnt_up_gene_list<-c("WNT3", "TCF7","TCF3","WNT9A","FZD9", "LEF1")
# graph settings
theme<-theme(panel.background = element_blank(),
panel.border=element_rect(fill=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background=element_blank(),
axis.text.x=element_text(colour="black"),
axis.text.y=element_text(colour="black"),
axis.ticks=element_line(colour="black"),
plot.margin=unit(c(1,1,1,1),"line"),
legend.title = element_text(face = "bold"),
legend.text = element_text(face = "bold"),
axis.title.x = element_text(color="black", size=18, face="bold"),
axis.title.y = element_text(color="black", size=18, face="bold"),
axis.text = element_text(size = 12, face = "bold"))
FC_df<-rbind(data.frame(genes =neuroendocrine_markers,logFC = DEGs_48h[neuroendocrine_markers,"logFC"], makers = rep("neuroendocrine_markers", n = length(neuroendocrine_markers))),
data.frame(genes =neural_genes, logFC = DEGs_48h[neural_genes,"logFC"], makers = rep("neural_genes", n = length(neural_genes))),
data.frame(genes =keratin_markers, logFC = DEGs_48h[keratin_markers, "logFC"], makers = rep("keratin_markers", n = length(keratin_markers))))
# Bar plot to show the FC of ER vs.GFP at 48 hrs
p<-ggplot(FC_df, aes(x=factor(genes, levels = genes), y=2^(logFC), fill=factor(makers, levels = c("neuroendocrine_markers", "neural_genes", "keratin_markers")))) +
geom_bar(stat="identity")+
guides(fill=guide_legend(title="Gene categories")) +
xlab("")+
ylab("Fold Change ER vs.GFP at 48 hrs") +
geom_hline(yintercept = 1, color = "red", linetype = "twodash") +
theme(legend.position="bottom")+
theme
p

# Display the statistic data from DEGs anaylsis for 48 hrs ER and GFP samples
DT::datatable(DEGs_48h[c(neuroendocrine_markers, neural_genes, keratin_markers),])
Plotting the WNT
genes expression dynamic pattern in ER samples (adjusted by GFP)
# dynamic expression pattern of GOIs
p_up<-ggplot(IMR90_ER_GFP_DELTA_df[IMR90_ER_GFP_DELTA_df$genes %in% wnt_up_gene_list, ],
aes(x=factor(time, levels = c(0,4,8,12,16,20,24,32,40,48)),
y=`ER-GFP`, fill = genes, colour = genes, group = genes)) +
geom_line(size =0.5)+
geom_point(size=1.5)+
geom_dl(aes(label = genes), method = list(dl.trans(x = x + 0.15), "last.points", cex = 0.7, face = "bold"))+
ylim(-3, 2)+
labs(x = 'time', y = "Relative gene expression level to control at the same time point")+
theme(plot.title = element_text(size = 12, face = "bold"),
legend.title=element_text(size=15),
legend.text=element_text(size=12))
p_up<-p_up+theme
p_down<-ggplot(IMR90_ER_GFP_DELTA_df[IMR90_ER_GFP_DELTA_df$genes %in% wnt_down_gene_list, ],
aes(x=factor(time, levels = c(0,4,8,12,16,20,24,32,40,48)),
y=`ER-GFP`, fill = genes, colour = genes, group = genes)) +
geom_line(size =0.5)+
geom_point(size=1.5)+
geom_dl(aes(label = genes), method = list(dl.trans(x = x + 0.15), "last.points", cex = 0.7, face = "bold"))+
ylim(-3, 2)+
labs(x = 'time', y = "Relative gene expression level to control at the same time point")+
theme(plot.title = element_text(size = 12, face = "bold"),
legend.title=element_text(size=15),
legend.text=element_text(size=12))
p_down<-p_down+theme
IMR90_WNT_genes_plot<-ggarrange(p_up, p_down,
ncol = 2,
nrow = 1
)
IMR90_WNT_genes_plot
